Molar volume of solid isotopic helium mixtures 
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Solid isotopic helium mixtures have been studied by path-integral Monte Carlo simulations in 
the isothermal-isobaric ensemble. This method allowed us to study the molar volume as a function 
of temperature, pressure, and isotopic composition. At 25 K and 0.2 GPa, the relative difference 
between molar volumes of isotopically-pure crystals of ^He and *He is found to be about 3%. This 
difference decreases under pressure, and for 12 GPa it is smaller than 1%. For isotopically-mixed 
crystals, a linear relation between lattice parameters and concentrations of helium isotopes is found, 
in agreement with Vegard's law. The virtual crystal approximation, valid for isotopic mixtures of 
heavier atoms, does not give reliable results for solid solutions of helium isotopes. 
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I. INTRODUCTION 

The lattice parameters of two chemically identical crys- 
tals with different isotopic composition are not equal, 
lighter isotopes giving rise to larger lattice parameters. 
This is due to the dependence of atomic vibrational am- 
plitudes upon atomic mass, along with the anharmonicity 
of the vibrations. This effect is most important at low 
temperatures, since the zero-point amplitude decreases 
with increasing atomic mass. At higher temperatures, 
the isotope effect on the crystal volume is less relevant, 
and disappears in the high-temperature (classical) limit 
at r > Qd (8d, Debye temperature), where vibrational 
amplitudes become independent of the atomic mass. At 
present, the isotopic effect in the lattice parameters of 
crystals can be measured with high precisioi*^. 

Other quantities, such as the vibrational energy, dis- 
play an isotope dependence at low T in a harmonic ap- 
proximation, due to the usual rescaling of the phonon 
frequencies with the isotopic mass M (w oc Af^^/^), but 
this dependence can show appreciable changes when an- 
harmonic effects are present. All these effects are ex- 
pected to be more important in the case of helium than 
for heavier atoms, fn fact, solid helium is in many re- 
spects an archetypal "quantum solid," where zero-point 
energy and associated anharmonic effects are apprecia- 
bly larger than in most known solids. This gives rise to 
peculiar properties, whose understanding has presented 
a challence for theories and modelling from a microscopic 
point of viewSi. 

Anharmonic effects in solids, and in solid helium in 
particular, have been studied theoretically for many years 
by using techniques such as quasiharmonic approxima- 
tions and self-consistent phonon theories^i^. An alter- 
native procedure is based on the Feynman path-integral 
formulation of statistical mechanics.^'-, that has turned 
out to be a convenient approach to study thermodynamic 
properties of solids at temperatures lower than their De- 
bye temperature Qd, where the quantum character of 
the atomic nuclei is relevant. In particular, Monte Carlo 
or molecular dynamics sampling can be applied to eval- 



uate finite-temperature path integrals, thus allowing one 
to carry out quantitative and non-perturbative studies of 
anharmonic effects in solids-. 

The path-integral Monte Carlo (PIMC) method has 
been employed to study several structural and thermody- 
namic properties of solid helium^'^'^ i^^'^"'^ , as well as heav- 
ier rare-gas solid a^^i^'^i-'^^'-'^^i"'^^'"'^'^ . For helium, in particu- 
lar, this procedure has predicted kinetic-energy values^ 
and Debye- Waller factor^i^ in good agreement with data 
derived from experimentsi^ii^. PIMC simulations have 
been also performed to study the isotopic shift in the 
helium melting pressur o^i^° . 

In most calculations of properties of crystals with iso- 
topically mixed composition, it is usually assumed that 
each atomic nuclei in the solid has a mass equal to the 
average mass. This kind of virtual- crystal approximation 
has been used in density-functional calculations, as well 
as in PIMC simulations— i^-^'^^i^^'^'^i^^ . In fact, in earlier 
simulations it was found that the results obtained by us- 
ing this approximation are indistinguishable from those 
derived from simulations in which actual isotopic mix- 
tures were considered. This seems to be true for atoms 
heavier than helium, and in particular for rare-gas solids 
including solid N©ii, but is not guaranteed to happen 
for solid helium, due to its low atomic mass and large 
anharmonicity. 

It is well kown that at temperatures lower than 1 K, a 
phase separation appears in solid '^He-^He mixtures, and 
the actual temperature at which this separation occurs 
depends on pressure and isotopic composition^LSl. This 
isotope segregation is due to the different molar volume 
of both isotopes, which in turn is caused by the different 
zero-point vibrational amplitudes of •^He and ^He^ii^^. 
In this paper, we consider mixtures of helium isotopes 
at higher temperatures, where '^He and *IIe form solid 
solutions for any isotopic composition. By varying the 
molar fraction of both isotopes, we analyze changes in 
the lattice parameter and kinetic energy, by using PIMC 
simulations. We employ the isothermal-isobaric (NPT) 
ensemble, which allows us to consider properties of these 
solid solutions along well-defined isobars. This Simula- 
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tion method permits to study properties of actual iso- 
topic mixtures, and compare them with those obtained 
for virtual crystals in which each atom has a mass equal 
to the average mass of the considered isotope mixture. 

The paper is organized as follows. In Sec. II, the com- 
putational method is described. In Sec. Ill we present the 
results, and Sec. IV includes a discussion and the conclu- 
sions. 



II. METHOD 

Equilibrium properties of solid helium in the face- 
centred cubic (fee) and hexagonal close-packed (hep) 
phases have been calculated by PIMC simulations in the 
NPT ensemble. Our simulations were performed on su- 
percells of the fee and hep unit cells, including 500 and 
432 atoms respectively. These supercell sizes are enough 
for convergence of the quantities studied here^S. For a 
given average isotopic mass Af , we randomly distribute 
■^He and "^He atoms in the appropriate proportions over 
the lattice sites of the simulation cell, and the atoms are 
kept fixed at their respective positions along the simula- 
tion (no diffusive positional changes). This is assumed to 
be valid for the temperatures studied here, much higher 
than those at which phase separation appears (T < 1 
j^-^j27_ Pqj. ga^(2]^ isotopic composition studied here, we 
have taken five different realizations of the isotope mix- 
ture. To analyze with more detail the dispersion in the 
lattice parameters obtained in the simulations, we took 
12 different samples for M — 3.25 amu, and 15 samples 
for M — 3.5 amu (see below). For comparison, we have 
also considered virtual crystals, in which every atom has 
a mass M, i.e., = M for i = 1, N. 

Helium atoms have been considered as quantum par- 
ticles interacting through an effective interatomic poten- 
tial, composed of a two-body and a three-body part. For 
the two-body interaction, we have taken the potential 
developed by Aziz et al?° (the so-called HFD-B3-FCI1 
potential). For the three-body part we have employed 
a Bruch-McGee-type potentia l^^'^^ , with the parameters 
given by Loubeyre^, but with parameter A in the attrac- 
tive exchange interaction rescaled by a factor 2/3, as in 
Ref J^. This interatomic potential has been found to de- 
scribe well the vibrational energy and equation-of-state of 
solid helium in the available range of experimental data, 
including pressures on the order of 50 GPa^^. 

The PIMC method relies on an isomorphism between 
the considered quantum system and a classical one, ob- 
tained by replacing each quantum particle by a cyclic 
chain of Ntt classical particles {Ntt- Trotter num- 
ber) , connected by harmonic springs with a temperature- 
dependent constant. This isomorphism appears because 
of a discretization of the density matrix along cyclic 
paths, that is usual in the path-integral formulation of 
statistical mechanics^'-^. Details on this computational 
method can be found elsewhere^i^^i^. 

Our simulations were based on the so-called "primi- 



tive" form of PIMG^^i^. We considered explicitly two- 
and three-body terms in the simulations, which did not 
allow us to employ effective forms for the density ma- 
trix, developed to appreciably simplify the calculation 
when only two-body terms are explicitly considered^*'. 
Quantum exchange effects between atomic nuclei were 
not taken into account, since they are negligible for solid 
helium at the temperatures and pressures studied here. 
(This is expected to be valid as long as there are no vacan- 
cies and T is higher than the exchange frequency ~ 10~^ 
K^.) For the energy we have used the "crude" estimator, 
as defined in Refsi^i^. 

Sampling of the configuration space has been carried 
out by the Metropohs method at pressures P < 12 GPa, 
and temperatures between 10 K and the melting tem- 
perature at each considered pressure. However, most of 
the simulations presented in this paper were carried out 
for fee He at 25 K and 0.3 GPa, conditions at which the 
isotopic effects studied here are clearly observable. Some 
simulations were also carried out at lower temperatures 
(see below). For given temperature and pressure, a typ- 
ical run consisted of 10* Monte Carlo steps for system 
equilibration, followed by 10^ steps for the calculation of 
ensemble average properties. To keep roughly constant 
the accuracy of the computed quantities at different tem- 
peratures, we have taken a Trotter number that scales as 
the inverse temperature 1/T. At a given T, the value 
of Nxr required to reach convergence of the results de- 
pends on the Debye temperature, higher requiring 
larger Nxr- Since vibrational frequencies increase as the 
applied pressure rises, Nxr has to be raised accordingly. 
For pressures on the order of 1 GPa, NttT = 2000 K is 
enough to reach convergence of the computed quantities. 
For pressures larger than 2 GPa, we have taken NttT — 
4000 K for ^He and 3000 K for "^He, as in earher workS^. 
Other technical details are the same as those employed 
in Refs^i^s^. 

In the isothermal-isobaric ensemble, the mean-square 
fluctuations in the volume V of the simulation cell are 
given by 



B ^ 



(1) 



where B — —V{dP/dV)T is the isothermal bulk 
modulus^. Hence, for a cubic crystal, the fluctuations 
in the lattice parameter a are: 



9L^aB 



(2) 



where is the number of unit cells in a simulation cell 
with side length La. From Eq. ^ one can see that the 
relative fluctuation in the volume of the simulation cell, 
ay/V, scales as L~^^'^. For *He at 25 K we found in 
PIMC simulations ay/V = 5.5 x 10"^ fo^. p ^ o.3 GPa, 
and 1.7x 10-3 for 12 GPa (for L = 5 and 500 atoms). For 
fee *He at 0.3 GPa this translates into da = 7.3 x 10"^ A. 
We are interested in an accuracy in a on the order of 10"* 
A, which means that at this temperature and pressure 
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FIG. 1: Convergence of the lattice parameter a of fee He as a 
function of the inverse Trotter number, A'^^^^ , as derived from 
PIMC simulations at T = 25 K and P = 0.3 GPa. Squares 
and circles correspond to "^He and *He, respectively. Dashed 
lines are guides to the eye. Error bars are smaller than the 
symbol size. 
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FIG. 2: Temperature dependence of the lattice parameter a 
of fee helium, derived from PIMC simulations at two differ- 
ent pressures: 0.3 and 0.6 GPa. Squares and circles indicate 
results for ^He and ''He, respectively. Error bars are smaller 
than the symbol size. Dashed lines are guides to the eye. 



one needs about 10^ independent data. In practice, the 
required number of Monte Carlo steps is expected to be 
larger, due to correlation between configurations along a 
Monte Carlo trajectory. In fact, we have checked that 
10^ Monte Carlo steps are enough to have an statistical 
uncertainty in a on the order of 10^^ A. 



III. RESULTS 

We have checked the convergence with the Trotter 
number of several quantities derived from our PIMC sim- 
ulations. In Fig. 1 we display the dependence of the lat- 
tice parameter a of fee '^He and ^He, as a function of 
the inverse Trotter number N^^ , for T = 25 K and P = 
0.3 GPa. The lattice parameter obtained in the simula- 
tions increases with Nxr, and converges to a finite value 
for large Nrr (limit N:^^ ^ in Fig. 1). The differ- 
ence 5a = as — a4 between lattice parameters of ^He and 
^He decreases as the Trotter number is lowered, and goes 
to zero in the classical limit (Nxr — 1), where this iso- 
topic effect disappears. The reliability of the interatomic 
potential employed here to predict lattice parameters of 
solid helium has been studied elsewhere^^. Here we will 
only comment that it gives good agreement with experi- 
mental results up to the melting temperature in the pres- 
sure range considered in this paper. Thus, for fee ^He at 
38.5 K and 0.493 GPa we find a = 3.9154 A vs 3.915(2) 
A derived from inelastic neutron scattering2£. 

Once checked its convergence with Nxr, in Fig. 2 we 
show the temperature dependence of the lattice param- 



eter a of fee '^He and ^He at two different pressures: 
0.3 GPa (open symbols) and 0.6 GPa (filled symbols). 
Squares correspond to '^He and circles to *He. For each 
pressure, results are displayed for temperatures at which 
the considered solids were found to be stable along the 
PIMC simulations. As expected, a is larger for '^He than 
for "^He, and the difference 5a = a^ — 04 is smaller for 
higher pressure. Also, for a given pressure, 5a decreases 
slowly as the temperature is raised (at higher T the solid 
becomes "more classical"). 

To quantify the change in crystal volume with isotopic 
mass, we calculate the ratio 5V/V4 = (V3 — T4)/V4. This 
ratio is shown in Fig. 3 as a function of pressure for fee 
(filled symbols) and hep He (open symbols) at 25 K. At P 
= 0.2 GPa, we find 5V/V4 = 0.030. This value is reduced 
by more than a factor of 3 at P = 12 GPa, where we ob- 
tain 5V/V4 = 7.7 X 10"'^. It is interesting to note that at 
P = 1 GPa, both fee and hep phases could be simulated 
at 25 K, remaining (meta)stable along the corresponding 
simulation runs. The ratio 5V/V4 obtained in both cases 
is shown in Fig. 3, and the results coincide within sta- 
tistical errors. Direct experimental measurements of this 
isotopic effect on the molar volume are scarce. Stewart^^ 
measured the molar volumes of both '^He and ''He at 4.2 
K and pressures up to 2 GPa. In particular, for P = 0.3 
and 0.6 GPa he found V4/V3 = 1.023. At these pressures, 
we found for the low-temperature limit in the PIMC sim- 
ulations V4/V3 = 1.026(1) and 1.023(1), respectively. 

The difference 5V is largest at small pressures and low 
temperatures, where quantum effects are most impor- 
tant. For a given solid, quantum effects on the crystal 
volume can be measured by the difference V—Vd between 
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FIG. 3: Isotopic effect on the crystal volume of solid he- 
lium, as obtained from PIMC simulations. Shown is the ratio 
SV/V4 — (V3 — V4) /V4 as a function of pressure at 25 K. Open 
and filled symbols correspond to hep and fee phases, respec- 
tively. Error bars are less than the symbol size. The dashed 
line is a guide to the eye. 
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FIG. 4: Kinetic energy of fee helium as a function of the mean 
isotopic mass M in isotopically mixed crystals. Symbols are 
results of PIMC simulations at 25 K and 0.3 GPa. The dashed 
line is a linear fit to the data points. 



the actual volume V and that obtained for a "classical" 
crystal of point particles, Vc\. This difference decreases 
for increasing atomic mass and temperatur o^'^'^^ . From 
our PIMC simulations at T = 25 K and a relatively low 
pressure of 0.3 GPa, we found an increase in the volume 
of solid ^He and ^He of 26% and 22% respectively, as 
compared to the "classical" crystal at zero temperature. 

PIMC simulations allow us to obtain the kinetic energy 
of the different atoms in the simulation cell at finite tem- 
peratures. For pure fee crystals of ^He and ''He, we find 
a kinetic energy Ek — 9.22 and 8.21±0.03 meV/atom, 
respectively (at 25 K and 0.3 GPa). This translates into 
a kinetic-energy ratio of 1.12, slightly lower than the ra- 
tio between zero-point energies expected in a harmonic 
approximation {Eq/Eq — ^^4/3 — 1.15). These energy 
values are similar to those obtained earlier from PIMC 
simulations in the NVT ensemble. In particular, our re- 
sults for "'He (giving a molar volume of 9.95 cm'^/mol) 
are in line with those obtained earlier in the NVT en- 
semble at temperatures close to 25 K and molar volumes 
around 10 cm^/mol'®, which in turn agree with experi- 
mental measurements''^. 

For isotopically-mixed crystals, one expects the kinetic 
energy Ek of the whole crystal to evolve smoothly as a 
function of the mean isotopic mass M . In Fig. 4 we show 
the kinetic energy of fee hehum vs M at 25 K and 0.3 
GPa, as derived from our PIMC simulations. Within 
the precision of our results, we observe a linear depen- 
dence of Ek vs M, as indicated by the dashed line. We 
have compared these results with those obtained from 
PIMC simulations for the virtual crystal with mass M, 



and found that differences between both sets of results 
are smaller than the statistical noise. A discussion on the 
relation between vibrational kinetic and potential ener- 
gies in solid helium was given elsewhere^, and will not 
be repeated here. 

Going back to the lattice parameter of fee He, in Fig. 5 
we display the dependence of a on the average isotopic 
mass, at 25 K and 0.3 GPa. Open squares correspond 
to simulations in which "^He and ^He atoms were ran- 
domly distributed over the crystal sites, according to the 
required average mass M . The results show a linear de- 
pendence of a on M . In Fig. 5 we also show results of 
PIMC simulations of fee He, where each atom has a fic- 
titious mass M (M^ = M for z =J_, TV). This virtual- 
crystal approximation yields for M different from 3 and 
4 amu a lattice parameter smaller than that obtained 
for a realistic distribution of isotopes over the simulation 
cell. For an isotopic mixture including 50% of each iso- 
tope, we find a = 4.06162 A vs Ave = 4.05950 A for the 
virtual-crystal approximation, i.e., a difference between 
both models a — a^c = 2.12 x 10""^ A. For this composi- 
tion we have considered 15 different isotope distributions, 
and found for the lattice parameter a standard deviation 
a = 2.5 X 10"'' A, a little smaller than the symbol size 
in Fig. 5. Thus, the difference a — a^c amounts to about 
8a. For 75% '^He, we took 12 realizations of the isotope 
mixture, and found a — a^c — 1-64 x lO^'' A. In this case, 
tj = 2.8 X 10~^ A, or a — a^c ~ 6(t. 
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FIG. 5: Lattice parameter a of fee helium as a function of 
the mean isotopic mass M. Data points were obtained from 
PIMC simulations at 25 K and 0.3 GPa, for isotopically mixed 
crystals containing ^He and ''He (squares), and for virtual 
crystals built up by atoms with the average isotopic mass 
(circles). Error bars are smaller than the symbol size. The 
dashed line is a least-square fit to the data points (squares). 
The dotted line was obtained from a quasi-harmonic approx- 
imation using Eq. (5). 



IV. DISCUSSION 



where C3 is the concentration of ^He. From the data 
presented in Fig. 5 we find an expansion coefficient 
(3 = 1.49 X lO^^'"^ cm'^/atom. According to our results, 
this linear relation between lattice parameter and isotope 
concentration holds for the whole concentration range, 
as shown in Fig. 5 (squares and dashed line). This is in 
agreement with earlier calculations for alloys, which indi- 
cate that for a difference between atomic radii less than 
~ 5% one does not expect appreciable departure from 
linearity^. In fact, ^He and ^He behave as atoms with 
a difference in effective atomic radii of about 1% (for T 
and P considered here). 

The virtual-crystal approximation has been employed 
earlier to study the dependence of lattice parameters on 
the average isotopic mass. This can be conveniently done 
by using a quasiharmonic approach, according to which 
the low-tcmperature lattice parameter a for mean iso- 
topic mass M can be approximated by^i 
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^«(q)7«(q) 



(4) 



Here, w„(q) are the frequencies of the nth mode 
in the crystal, B is the bulk modulus, doo is the 
zero-temperature lattice constant in the limit of in- 
finite atomic mass (classical limit), and 7n(q) — 
—d\nu!n{ci) /dlnV is the Griineisen parameter of mode 
n, q. Assuming a mass dependence of the frequencies 

Wn(q) 

tice parameter with isotopic mass 



1/2 

M , one finds for the relative change in lat- 



Path-integral Monte Carlo simulations have been 
found to be well suited to study finite-temperature anhar- 
monic effects on structural and thermodynamic proper- 
ties of crystalline solids. These effects are particularly 
important for solid helium, where isotopic effects are 
relevant, as manifested in differences in the molar vol- 
ume and vibrational energies of solid '^He and ^He. The 
PIMC method enables us to study phonon-related prop- 
erties without the assumptions of quasiharmonic or self- 
consistent phonon approximations, and to study anhar- 
monic effects in solids in a non-perturbative way. Thus, 
for a given reliable interatomic potential, this method 
yields in principle "exact" values for measurable proper- 
ties of many-body quantum problems, with an accuracy 
limited by the imaginary-time step (Trotter number) and 
the statistical uncertainty of the Monte Carlo sampling. 

Our results for the lattice parameter of solid helium 
as a function of the average isotopic mass (see Fig. 5) 
can be understood in terms of the theory of alloys. In 
fact, ^He and ^He behave in this respect as atoms with 
different atomic radii, as a consequence of the different 
vibrational amplitudes of both isotopes. If we consider 
a ''He crystal with '^He impurities, the lattice expansion 
Aa due to these impurities is given by Vegard's law^^iii: 



AM 



(5) 



Aa 
04 



PCs 



(3) 



where A is constant for a given pressure. By applying this 
equation to 03 and 04, we can find A, and then the de- 
pendence of a on M shown in Fig. 5 as a dotted line. This 
line coincides within error bars with the results derived 
from the PIMC simulations for helium with Mi = M for 
all i (virtual crystal) . Note that the low-temperature ex- 
pression for the lattice parameter given in Eq. (4) is a 
good approximation for the conditions (P, T) considered 
here. In fact, for a volume of ~ 10 cm^/mole, a tem- 
perature of 25 K can be considered a "low temperature" 
when compared with 0_d ^ 120 K— . 

It is worth commenting on the difference found be- 
tween this quasiharmonic approximation and the actual 
isotope distribution over the crystal. For the virtual crys- 
tal we found in both PIMC simulations and in the quasi- 
harmonic approach, a non-linear dependence of the lat- 
tice parameter a on average isotopic mass. This contrasts 
with the linear dependence derived from PIMC simula- 
tions for actual distributions of isotopes. (If a depar- 
ture from linearity appears in this case, it will be smaller 
than the precision of our results.) This observation is 
important for the helium solid solutions considered here, 
since in most known solids both approaches give the same 
resultsi^i^ii^. From the results displayed in Fig. 5, we 
observe for M = 3.5 amu a difference a — Qvc = 2.1 x 10""^ 
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A, which amounts to about 0.05% of the lattice parame- 
ter. This relative difference is much larger than the un- 
certainty in structural parameters currently derived from 
diffraction methods^. 

In the last few years, several authors have indicated 
that pressure causes a decrease in anharmonicity^^^i^, 
in agreement with earlier observations that the accuracy 
of quasiharmonic approximations increases as pressure is 
raised and the density of the solid increases^i^. This 
is also the origin of the decrease in the isotopic effects 
studied here, as pressure is raised. Since these isotopic 
effects are caused by anharmonicity of the interatomic in- 
teractions, an effective decrease in anharmonicity causes 
a reduction in the difference Sa = 03 — 04 or, equivalently, 
in the ratio 5V/V4 shown in Fig. 3. 

In summary, we have carried out PIMC simulations 
of solid solutions of helium isotopes in the isothermal- 
isobaric ensemble. Our results indicate that Vegard's law 



is fulfilled in the whole composition range, i.e., the crys- 
tal volume changes linearly with isotopic composition. 
This volume change decreases appreciably as pressure is 
raised, but it is still clearly observable at pressures of the 
order of 10 GPa. Approximations such as virtual crystals 
with average atomic mass are not valid for solid isotopic 
helium mixtures. 
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